%This file generates a case of a 20x20 Cov-matrix with up to 10 sub blocks.
%sampled sequences with different levels of noice are generated and used 
%as estimated subblocks. The performance is measured for different number 
%of blocks and different noise levels. Sample length is fixed to 5000.

clear, close all
NumOfStocks = 20;

load data_10_blocks.mat

TestBlocks=[3 4 5 6 7 8 9 10];
%TestBlocks=[3 8];
%TestBlocks=[3:10];

error=zeros(length(TestBlocks),1);

for iterBlocks=1:length(TestBlocks)

    NumOfBlocks=TestBlocks(iterBlocks);
    
    %Define weights (dummy)
    w=ones(1,NumOfBlocks);
    
    
    %Generate the optimal estimate
     cvx_begin
        cvx_quiet(true); 
        variable Sigma_hat(NumOfStocks, NumOfStocks) symmetric;

        %Define objective function
        f = 0;
        s_max = 0;
        for q=1:NumOfBlocks
            f = f + w(1,q)*norm(Sigma_hat(Combinations{q}, Combinations{q})-CovBlocks{q},'fro');           
        end
        minimize (f)

        subject to

        Sigma_hat == semidefinite(NumOfStocks)
        
        
    cvx_end
    
    %Compare estimate to original matrix
    error(iterBlocks)=norm(Sigma-Sigma_hat,'fro');
    DiffMatrixRel(:,:,iterBlocks)=abs((Sigma_hat-Sigma)./Sigma);
    DiffMatrixAbs(:,:,iterBlocks)=abs((Sigma_hat-Sigma));
end

% figure(1)
% for i=1:length(TestBlocks)
% subplot(1,2,i)
% clims=[0 0.12*max(max(max(abs(DiffMatrixRel(:,:,:)))))];
% imagesc(DiffMatrixRel(:,:,i),clims);
% colormap(gray)
% colorbar
% title(['Relative estimation error with ', num2str(TestBlocks(i)), ' sub blocks']) 
% 
% end

figure(3)
plot(TestBlocks,error,'-k','linewidth',3)
title('Estimation error')
xlabel('Number of subblocks used')
ylabel('||S-Shat||_f')
set(get(gca,'title'),'fontsize',14)
set(get(gca,'xlabel'),'fontsize',14)
set(get(gca,'ylabel'),'fontsize',14)


% figure(2)
% for i=1:length(TestBlocks)
% subplot(2,2,i)
% clims=[0 0.5*max(max(max(abs(DiffMatrixAbs(:,:,:)))))];
% imagesc(DiffMatrixAbs(:,:,i),clims);
% colormap(gray)
% colorbar
% title(['Absolute estimation error with ', num2str(TestBlocks(i)), ' sub blocks']) 
% set(get(gca,'title'),'fontsize',14)
% end

% %Get correlationcoefficients
%     for i=1:NumOfStocks
%         for j=1:NumOfStocks
%             SigmaCoeff(i,j)=Sigma(i,j)/(sqrt(Sigma(i,i)*Sigma(j,j)));
%             Sigma_hatCoeff(i,j)=Sigma_hat(i,j)/(sqrt(Sigma_hat(i,i)*Sigma_hat(j,j)));
%         end
%     end
%     
% figure(2)
% for i=1:length(TestBlocks)
% subplot(2,2,i)
% clims=[0 0.5*max(max(max(abs(DiffCoeffMatrix(:,:,:)))))]
% imagesc(DiffCoeffMatrix(:,:,i),clims);
% colorbar
% end
